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Abstract 

In most high energy cosmic ray surface arrays, the primary energy is currently 
determined from the value of the lateral distribution function at a fixed distance 
from the shower core, r$. The value of tq is mainly related to the geometry of the 
array and is, therefore, considered as fixed independently of the shower energy or 
direction. We argue, however, that the dependence of ro on energy and zenith angle 
is not negligible. Therefore, in the present work we propose a new characteristic 
distance, which we call r op t, specifically determined for each individual shower, with 
the objective of optimizing the energy reconstruction. This parameter may not only 
improve the energy determination, but also allow a more reliable reconstruction of 
the shape and position of rapidly varying spectral features. We show that the use 
of a specific r op t determined on a shower-to-shower basis, instead of using a fixed 
characteristic value, is of particular benefit in dealing with the energy reconstruction 
of events with saturated detectors, which are in general a large fraction of all the 
events detected by an array as energy increases. Furthermore, the r op t approach has 
the additional advantage of applying the same unified treatment for all detected 
events, regardless of whether they have saturated detectors or not. 
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1 Introduction 



The energy spectrum of cosmic rays (CR) is observed from energies below 
1 GeV up to more than 10 20 eV. CR hit the top of the atmosphere at a 
rate > 10 3 per square meter per second. For energies below ~ 10 2 TeV, CR 
can be detected by direct measurement from high altitude balloons or satel- 
lites, albeit with rapidly decreasing statistics since the CR spectrum is very 
steep, ~ E~ 2,7 . For E > 10 2 TeV, the point of maximum development of 
the cascades of daughter particles initiated by collisions with atmospheric nu- 
clei starts reaching the ground at high altitudes. From that point onward, 
the properties of primary cosmic ray can be determined indirectly from the 
measurement of extensive air showers (EAS). 

Two different experimental approaches have been traditionally used to study 
the highest energy EAS. The first one consists on the inference of lateral 
distribution of particles from the discrete sampling of the shower front at 
ground level by a surface array of detectors. Scintillators (e.g., Volcano Ranch, 
AGASA, KASCADE) and water Cherenkov detectors (e.g., Haverah Park) 
have been mainly used for this purpose. The second technique consists on the 
reconstruction of the longitudinal profile of the shower from the fluorescence 
light produced by atmospheric Nitrogen as it is excited by charged EAS par- 
ticles along the atmosphere (e.g., Fly's Eye, HiRes). The latter is considered 
to be close to a calorimetric measurement of the primary CR particle energy. 
Extensive reviews on CR theory and experiments can be found in [1,2,3]. 

A special case from the experimental point of view is the Pierre Auger Ob- 
servatory [4] which pioneers the simultaneous use of both techniques, water 
Cherenkov detectors and fluorescence telescopes. For these hybrid events, sys- 
tematic errors in their energy estimate are greatly reduced. Unfortunately, 
fluorescence can only be observed during dark nights and, consequently, this 
technique can only be applied to ~ 13% of the incoming events. Therefore, 
even if the hybrid technique, when simultaneously available, allows for an in- 
dependent calibration of a ground detector, high statistics must come from 
surface arrays with a duty cycle close to 100%. 

In this work we are interested on the highest energies, E > 10 17 eV, at which 
point extragalactic flux is likely to penetrate the Galaxy and start dominating 
the CR flux. At these energies at least three very important spectral features 
are located: the second knee, the ankle and the GZK complex (bump plus steep 
flux suppression) [5]. The determination of their exact position and shape is 
a fundamental experimental problem. Therefore, the most relevant parameter 
in this spectral region is, arguably, the primary energy of the impinging CR. 

The procedure to determine the primary energy in sufrace arrays is a two step 
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process. First, the lateral distribution function (LDF), i.e. the shower particle 
density or signal versus distance to the shower axis, is fitted assuming a known 
functional form. This fit suffers from uncertainties related to the statistical 
shower fluctuations, the uncertainties in core location and the ignorance of 
the exact form of the LDF. The normalization constant of the LDF of an 
extensive air shower is a monotonous (almost linear) increasing function of 
the energy of the primary cosmic ray. Therefore, Hillas [6] proposed to use 
the interpolated signal at some fixed, characteristic distance from the shower 
core, S(r ), at which fluctuations in the LDF are minimal. The uncertainty 
due to the lack of knowledge of the LDF is also minimized by this procedure 
[7]. The use of the signal interpolated at r , S(r ) is widely used as energy 
estimator by surface detector arrays. AGASA [8,9], Yakutsk [10] and Haverah 
Park [11], for example, choose r = 600 m, while Auger uses 1000 m due to 
its larger array spacing [12]. The characteristic distance r is mainly, although 
not completely, determined by the geometry of the array. Thus, the same 
value of ro is used to estimate the energy for all the showers, regardless of 
primary energy or incoming direction. In the second step, there are at least 
two possible approaches to calibrate S(r ) as a function of primary energy: 
either via Monte Carlo simulations or, as in the case of Auger, by using the 
almost calorimetric measurement obtained from the fluorescence observation 
of high quality hybrid events [12] . 

As an alternative, but motivated by Hillas' original idea [6], in the present 
work we focus in the shower-to-shower determination of an optimal distance 
to the core, which we name hereafter r opt , at which the interpolation of the 
reconstructed signal is the best energy estimator for each individual shower, 
regardless of whether this point is actually the one that minimizes shower to 
shower LDF fluctuations. 

We perform a detailed study of r opt as function of array spacing, primary 
energy and the zenith angle of the incoming cosmic ray and demonstrate that, 
although array geometry is an important underlying factor, the dependence of 
r opt on the remaining parameters is not negligible. We study the bias associated 
with both techniques, r and r opt , and show that, if the dynamical range of 
the detector covers a wide interval of energies, it is much safer to estimate 
an r opt for the energy reconstruction of each individual event than to fix a 
single ro for the whole data set. In fact, not only the bias as a function of 
energy can be kept negligible over at least 2.5 decades in energy, but also the 
error distribution functions are much better behaved, i,e, without appreciable 
kurtosis or skewness and very much Gaussian in the mentioned energy range. 
The latter has a potential impact in the reconstruction accuracy of the energy 
spectrum. We demonstrate this by applying a fixed ro as well as a shower-to- 
shower r opt , to a simplified version of the actual energy spectrum between ~ 1 
and ~ 100 EeV. 
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A further advantage of the r opt approach is the straightforward treatment of 
events with saturated detectors. The problem of saturation is very common 
in all surface experiments, specially when dealing with high energy vertical 
showers. In fact, at the highest energies inside the designed dynamical range 
of any experiment, usually events with saturated detectors can account for a 
large, if not dominant, fraction of all the observed events. Different strategies 
have been used to deal with them. In some cases saturated detectors are 
directly discarded from the LDF fit, while in others the saturation value is 
used as a lower limit to the true signal during the fitting procedure. The Auger 
Collaboration is developing at present special, more sophisticated algorithms 
to estimate the signal of a saturated detector [13] in order to more properly 
account for them in the LDF fit. We show here that it is actually not possible to 
define a single characteristic r distance for both kinds of events. In fact, even 
if well defined medians values of r opt for events with and without saturated 
detectors do exist, the dispersions around the median at any energy are so 
large that both sets cannot be clearly differentiated as to use, for example, 
just two fixed distances instead of a single one. Nevertheless, using a shower- 
to-shower r opt distance, the inferred energy is unbiased for events with and 
without saturated detectors alike. This reconstruction strategy allows for an 
homogeneous treatment of the data set regardless of the increasing number 
events with saturated detectors when the energy increases. 

In a recent work [14], Newton and co-workers also estimated an optimal 
shower-to-shower distance, but used a different algorithm and with a some- 
what different scope. They were mainly concerned with demonstrating the 
existence of a single distance for any given shower at which fluctuations in the 
LDF are minimum. By assuming that such fluctuations can be well described 
by the fluctuations of just one parameter, the slope of the LDF, externally 
fed into their procedure, and using a combination of simulations and semi- 
analytical analysis, they claim that, regardless of the functional form of the 
LDF considered, there exist a convergence point of the LDFs, at a charac- 
teristic distance they call optimal, where shower-to-shower fluctuations are 
minimal. Their results, combined together for a mix of energies drawn from 
a flat spectrum, seem to support their claim and lead them to the conclusion 
that a single fixed distance, depending only on the geometry and spacing of a 
given array, would be a good choice for the energy determination in the whole 
energy range of the experiment. Furthermore, it is not clear from their study 
how to deal with the events with saturated detectors in the later scenario. 

Alternatively, in the present work we do not constrain the parameters of the 
LDF, which are an output of the fit to the simulated data. We introduce 
instead reliable error estimations for the reconstruction of the core position, 
as calculated by [15] for arrays of varying spacing as a function of energy. 
Furthermore, our final scope is the determination of energy all along the dy- 
namical range of an experiment, and not the study of the manifestation of 
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signal fluctuations in the LDF. Therefore, we study in detail the dependence 
of r opt and of its distribution function as a function of energy, zenith angle and 
array spacing. This study is performed for events with and without saturated 
detectors. We also give a comparative description of error and biases for the 
fix distance and the r opt distance approaches in that parameter space. In the 
same line, we further extend our analysis to the reconstruction of a simulated 
energy spectrum of known shape, and show what the potential effects are of 
using each technique. 

The paper is organized as follows. Section 2 describes our general algorithm. 
Two different detector arrays are considered, scintillators and water Cherenkov 
tanks. In Section 3, in order to study the r opt dependencies with array spacing 
and the energy and incoming direction of primary cosmic ray, water Cherenkov 
(Auger-like) stations have been used. In Section 4 we deal with the issue of 
energy determination. In that analysis, (AGASA-like) scintillators are con- 
sidered. A general discussion and conclusions are given in Section 5. While 
different detectors are used in Sections 3 and 4, the algorithm to find r opt is 
the same for both and the results and conclusions of the paper are not affected 
by the array under consideration. 



2 Algorithm 

The basic idea of our algorithm is to estimate the optimum distance to the 
core at which to determine the energy of a shower under the most realistic 
possible conditions. 

We assume a certain analytical LDF as the intrinsic average lateral distribu- 
tion of particles inside the shower front as a function of the distance to the 
core. For the chosen energy and geometry of the event (azimuth, zenith and 
true core position), this function is used to estimate the average LDF value 
at the actual position of each detector. Afterward, a signal is calculated using 
the previous average value as the mean of a Poissonian distribution. If the 
calculated signal falls between a minimum threshold and an upper limit corre- 
sponding to a saturation condition, it is assigned to the detector. In the case 
of saturation, the event is kept, but with a flag indicating this fact, although 
the saturated detectors are not used in the subsequent analysis, i.e. in fitting 
the LDF. Once a set of triggered detectors participating in the event has been 
defined, the reconstructed LDF is emulated by fitting an experimental LDF, 
which depends on the detector array under consideration, and is not neces- 
sarily the real LDF used in the first step to generate the event. The LDF 
fit requires an estimate of the core position. Such estimate is an important 
component of the analysis of the event and comes, in practice, from a global 
reconstruction procedure which implies an energy dependent error in the in- 
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ferred position of the core. In our algorithm, we simulate this error by shifting 
the reconstructed core position according to its experimentally determined 
Gaussian distribution function. For each shifted core position an independent 
LDF fit is performed. We define the optimum distance to the core r opt as the 
interpolated distance at which the dispersion between the several LDF fits 
is minimal. We argue that the interpolated signal at this point is the opti- 
mum estimator of the energy of a real event and constitutes the operational 
definition of our parameter r opt [16]. 

We use the following numerical approach to simulate EAS detection in a sur- 
face array. The array is a set of equally spaced detectors, located at the vertices 
of an infinite grid of triangular elementary cells with variable spacing. The in- 
put parameters of an event are its energy, azimuthal and zenith angles and core 
position. The identity of the primary particle is not taken into account since 
differences in composition produce only a small effect in the error distribution 
function of the reconstructed core position [15] which, in turn, when combined 
with the use of an experimental LDF maps into a negligible variation in both 
r and r opt . 

Whenever we simulate a water Cherenkov detector, we assume that the true 
lateral distribution of the signal is best represented by a Nishimura-Kamata- 
Greisen (NKG) function. This functional form was first obtained in an analyt- 
ical study of the lateral development of electromagnetic showers in [17], and 
later extended to the hadronic initiated showers because the electromagnetic 
particles represent around the 90% of the total particles of the shower. The 
NKG selected is normalized at 1 km in the same way as the reported by Auger 
in [4]: 



7 co p0.95 o/3(0) 

S(r, E, 6) = , '■ b4h 1 x x (1 + r )~M (1) 

ll + 11.8[sec(0) - l] 2 



where r is the distance to the shower axis expressed in km, E is the energy of 
the primary in EeV, 9 is the zenith angle and j3(9) = 3.1 — 0.7sec(8). The signal 
in eq. 1 is expressed in vertical equivalent muons ( VEM) , which correspond to 
the signal deposited by one vertical muon in an Auger water Cherenkov tank. 

We use eq. 1 as the real LDF to simulate any given incoming event. The 
measured signal at each station is obtained with a Poissonian probability 
distribution function whose mean is given by eq. 1, the "true" LDF. The 
trigger condition is set to S(r) = 3.0 VEM. The saturation value is fixed at 
S(0.2 km, 1 EeV, 0°). These values are compatible with the equivalent Auger 
parameters. 

The uncertainty in core determination depends on the array geometry and 
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primary energy and it has been estimated for a variety of cases in reference 
[15]. We simulate the reconstruction uncertainty of the core using a Gaussian 
distribution function centered at the position of the real core, with standard 
deviation given by [15] as a function of the energy of the shower for the array 
spacing under consideration. 

For any shower, the following procedure is used to obtain the optimum distance 
r opt . Throughout the procedure, we try to mimic, as far as possible, the actual 
reconstruction procedure. As explained earlier, several fits to the LDF are 
performed for any event, each one with its own estimated core position. Since 
the exact functional form of the LDF function is not crucial [7] we use a generic 
LDF parametrization to fit the signals of the triggered stations: 

logS(r) = a ir - a * +a 3 (2) 

The uncertainty in the core position used for each one of the LDF fits corre- 
sponding to a given event, is accounted for by randomly shifting that point 50 
times with the same Gaussian probability distribution function referred above 
centered at the position of the reconstructed core. 

For each new core position, the LDF fit is performed using eq. 2. The slope 
and the normalization constant of each LDF are determined from the fitting 
procedure. The r opt value is defined as the point at which the dispersion among 
the interpolated signals over the several LDFs goes through a minimum. 

Therefore, the implementation of the algorithm requires a two step process: 
first, a global fit to the LDF is performed in order to get an estimate of 
the reconstructed core position and, second, the reconstructed core itself is 
fluctuated and r opt obtained using the procedure explained in the previous 
paragraph. 

When simulating scintillators as those of the AGASA experiment, we follow 
exactly the same procedure as before, but we use the following LDF [18] 
instead of eq. 1: 

p(r, E, 6) = 49.676 x 7.55" (e) - L2 x f s (0) x E 1/im x (3) 
\r M J V r M J \ l 

where p is given in m -2 , distances are in m, tm = 91.6 m is the Moliere radius 
at AGASA altitude, r}{6) = 3.84-2.15(sec(6 l ) - 1) and f s (6) is the attenuation 
curve: 



7 



fs(0) 



e Xp [--(sec9 



1) 



A 2 



{secO-lf 



(4) 



where X Q = 920 g/cm 2 , Ai = 500 g/cm 2 and A 2 = 594 g/cm 2 for showers with 
< 45°. The signal is fluctuated, as always, with a Poissonian distribution. 

In this case, the trigger condition is selected in such a way that the signal 
is not dominated by fluctuations. In particular, we use a vale of p such that 
fluctuations account at most by 50% of the signal. The saturation value is 
p(0.2 km ,1 EeV ,0°). 

In the AGASA case, we use as shower generator eq. 3, and perform the sub- 
sequent fitting procedure using an LDF with the functional form "observed" 
by AGASA: 



which is formally equivalent to eq. 3 for r >> tm [8]. 



3 r opt dependence on array spacing, energy and zenith angle. 

We consider in this section water Cherenkov detectors with separations of 433, 
750, 866 and 1500 m, as well as primary energies varying from 10 17 to 10 19,5 
eV. We use eq. 1 to generate the signals and eq. 2 to fit the LDF. 

In all cases we consider a uniform distribution in azimuth and zenith angles 
as explained in each figure. Shower cores are uniformly distributed inside an 
elementary cell of the array. 

Events with and without saturated detectors lead frequently to systematically 
different behaviors regarding the relationship between r and r opt under dis- 
cussion and, in principle, should be treated differently during data processing. 
Thus, in what follows, we will analyze them separately whenever appropriate. 

Figure 1 shows the dependence of r opt with energy without discriminating 
whether events have saturated detectors (labeled as All). Showers are injected 
with zenith angles 9 = 30° and 60°. It can be seen that r opt is a monotonous 
increasing function of energy due to the triggering of stations progressively 
further away from the core as energy increases. In Figure 2 the same behavior 
is shown separately for events with and without saturation. 

Since events with and without saturated detectors may, and indeed do, behave 



logp(r) = ai — a 2 log(r/rjvf) — 0.61og(l + (r/1000m) 2 ) 



(5) 
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Fig. 1. r op t vs. energy for different array spacing and zenith angle. The error bars 
represent the 68% and 95% C.L. The labell All means that events with and without 
saturated detectors are both included. 

in different ways, Figure 2 shows r opt results for both separately. All the previ- 
ous array spacings are considered but only at one zenith angle, 9 = 30°. r opt is 
an increasing function of energy for both sets of events. It can be seen that r opt 
is larger for events with saturated detectors at lower energies (r^l > r™°™~ sat ) 
but that, at higher energies, Top/ 1-50 * rapidly grows towards r^®*. The transition 
energy region is narrow (AlogE ~ 0.2) and shifts upwards in energy as the 
array spacing grows. The symmetry of the triangular array with respect to a 
shower core located at the center of an elementary triangle (the least likely 
configuration to saturate at any given energy), manifest itself in ring-like ar- 
rangements of triggered stations. The appearance of a third ring of triggered 
stations is responsible for the rapid grow of r™r sat over a limited energy 
interval as shower energy grows. 

Furthermore, low energy events with saturation have their cores very near the 
saturated stations. Therefore, the first triggered stations that do not saturate 
are clustered at the same distance from the core, which is roughly the array 
separation distance. Therefore, it is at the array separation distance that the 
dispersion among the several fits to the LDF is minimum. At higher energies, 
however, the next ring of the array enters into the set of triggered detectors of 
the event and, naturally, r opt increases. In Figure 2 it can be seen that r opt is 
almost constant and very near to the array separation at the lower energies, 
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Fig. 2. r op t vs. energy. Events with and without saturated detectors are shown 
separately. Zenith angle is 9 = 30°. The error bars represent the 68% and 95% C.L. 

and that there is a threshold energy, which depends on the array separation, 
from which r opt increases steadily with energy. 

Figure 3 shows the dependence of r opt with zenith angle for the same array 
spacings and three different input energies: log(E/eV) = 18.5,19.0 and 19.5. 
Both, events with and without saturated detectors are included. It can be 
seen that r opt is almost independent of zenith angle for 9 < 30° for any array 
spacing. However, as the zenith angle increases beyond 30°, r opt decreases with 
9, independently of array spacing and energy. The same effect is observed in 
both sets of events, those with saturation (Figure 4. a) and without it (Figure 
4.b). This result comes from the fact that, for inclined showers, the array 
spacing projected onto the shower front shrinks with zenith angle and r opt 
naturally follows this behavior. 

From the previous results, it is clear that r opt is in general a function of energy 
and zenith angle for inclined showers. Furthermore, in Figures 1, 2, 3 and 4, 
the error bars indicate the 68% and 95% confidence levels (CL) and the central 
points correspond to the median value of r opt . It can be seen that, in all cases, 
even if the behavior of the median curves is rather smooth, the CL are large 
and, therefore, considerable fluctuations are expected. Additionally, due to 
the large relative fluctuations of the signals from detectors located at large 
distances from the shower axis, the error distributions are skewed in general 
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Fig. 3. r op t vs. zenith angle for different array spacing and energies. The error bars 
represent the 68% and 95% C.L. Events with and without saturated detectors are 
both included. 



"1.8 
a. ~ 

J 16; 
1.4 : 
1.2: 
1- 
0.8. 
0.6 
0.4 
0.2" 
0- 



J 8 ' 



_ . logE = (9.5 
• logE =19.0 



10 20 30 40 50 



+ logE* 


19.5 


+ logE = 


19.0 


. logE' 


18.5 




e [deg] 



B [deg] 



Fig. 4. r op t vs. zenith angle for 1500 m separation array. The error bars represent 
the 68% and 95% C.L. (a): Events with saturated detectors, (b): Events without 
saturated saturated detectors.. 

towards larger values from the median of r opt . These points argues strongly 
in favor of an r opt determined specifically for each shower since, using a fixed 
characteristic value, ro, could compromise the estimation of primary energy. 
This possibility is analyzed in the following section. 

As it was mentioned in the introduction, although similar in character, the 
work in reference [14] is rather different in algorithmic approach and scope. 
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Therefore, a comparison between results in both works is not straightforward. 
Nevertheless, Figures 5 in [14] can be used to some extent to crosscheck our 
results. Figure 5 bottom-right in [14] shows their r opt as a function of energy. 
Despite the fact that there are indications of border effects at low energies in 
their calculation and that different zenith angle events are binned together, 
the results are similar to those in our Figure 2.d. Figure 2 shows r opt for 
events with and without saturated detectors in the energy interval between 
~ 10 and 30 EeV, for 433 (a), 750 (b), 866 (c) and 1500 m (d) spacing. It can 
be seen that, at 433, 750 and 866 m spacing r opt is more or less independent 
of energy at lower energies but eventually increases steadily above a certain 
energy. This effect is also expected at a separation of 1500 m for energies 
beyond those presently plotted in Figure 2.d. Reference [14], on the other 
hand, shows results only for arrays at 1500 m separation, where the same 
trend seems to be suggested for events with saturated detectors (see Figure 5 
bottom-right of [14]). Remarkably, although their analysis extends up to 100 
EeV, the same trend is not seen for events without saturation. The latter, 
however, may be due to the fact that in reference [14] showers with all zenith 
angles are mixed together which, at high energies, implies that their sample 
must be highly biased to very inclined events (otherwise they would present 
saturation), masking the effect. In fact, it can be seen from our Figure 4.b 
that, for events without saturation, r opt does decrease at any energy for larger 
zenith angles. 

Again, in Figure 5 bottom- left of [14], and despite the fact that the authors 
claim only a slight dependence of r opt with zenith angle, we obtain a very 
similar result for events with saturation in Figure 4. a with r opt decreasing 
markedly with increasing zenith angle. There is no agreement, however, for 
events without saturation, where they obtain an r opt that increases with zenith 
angle, while our results (see, Figure 4.b) shows an r opt that at low energies 
decreases as a function of zenith angle, but tends to a constant value as the 
energy increases. Part of the difference between both results may be due to the 
fact that in [14] energies randomly selected from a flat spectrum have been 
binned together. The latter, however, cannot account for their unexpected 
raise with zenith angle. 



4 Influence of r op ( on reconstructed energy 

In this Section we analyze the effect of adopting a fix characteristic distance, 
r , instead of a shower-specific value, r opt , for the determination of shower 
energy and energy spectrum. We simulate a detector similar to AGASA (see 
Section 2), i.e., a separation of 1 km between stations and use eq. 3 and eq. 5 
in order to generate signals and fit the "observed" LDF respectively. For each 
event r opt is estimated using the procedure explained in Section 2 while eq. 3 
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is used in order to estimate the energy for both ro = 600 m, as AGASA did, 
and r opt . 

Two different input spectra are used. A spectrum with one thousand events 
per energy bin (Alog(E') = 0.1) from 10 17 ' 8 to 10 20,7 eV, is used in order to 
study the functional form of the energy error distributions and the energy re- 
construction bias (Sections 4.1 and 4.2). The energy reconstruction of events 
with saturated detectors is also analyzed. Second, in Section 4.3, a more struc- 
tured spectrum extending from 10 17,7 to 10 20 ' 5 eV, which possesses an ankle, 
a GZK-cut-off, and is exposure-limited at low energy, is used to assess the 
impact of both techniques in a more realistic situation. For every event, the 
angular distribution is extracted randomly from an isotropic distribution with 
a maximum zenith angle of 45°, as in the case of the AGASA experiment. The 
azimuthal angles are selected from a uniform distribution between 0° and 360° 
and the core location is randomly located inside an elementary cell. 

It must be noted that the results of this Section do not directly apply to the 
spectrum inferred from surface arrays that relay on the use of hybrid events 
for the energy calibration. The main reason is that the error in core location 
for hybrid events is much smaller than for pure surface events. In the case of 
Auger, for example, the error in hybrid core determination is only around 35 
m [19] while for the majority of pure SD events it is estimated to be around 
100 m [20]. Therefore, r opt for hybrid events is very much constrained. In addi- 
tion, hybrid experiments do not directly relate the signal measured at r opt with 
the primary energy. Instead, they use a calibration with the energy obtained 
by the fluorescence technique. Finally, the most important uncertainties in 
the primary energy determination in hybrid experiments come from the flu- 
orescence uncertainties not from the parameter size determination as will be 
discussed later in detail. 



4-1 Error distribution function in energy reconstruction 

We calculate the distribution function of the errors in reconstructed energy, 
i.e. the difference between the reconstructed and the real energy, as a function 
of injected energy for both techniques, r and r opt . 

Figure 5 shows, for both r (a) and r opt (b), the 68% and 95% CL for the 
right and left sides with respect to the median of the energy error distribu- 
tion. It can be clearly seen that the error distribution functions originated by 
using r opt are much more compact and symmetrical than the corresponding 
distributions for ro- The effect is more notable for lower energies where the 
distribution function for characteristic distance determination is particularly 
wide and skewed. Although these figures are drawn for the 1000 m separation 
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LogE(ieal]/eV LogE(ieal]/eV 



Fig. 5. 68% and 95% CL over the median, from both its lower and higher energy 
sides, for the energy error distribution functions determined using either tq (a) or 
r opt (b) methods. See text for more details. 

array, the results apply qualitatively for any of the other spacings considered 
in previous sections. 

Arguably, it is desirable that the errors in energy reconstruction have a Gaus- 
sian distribution. Gaussian errors, for example, are easier to handle and under- 
stand when applying deconvolution techniques in the spectrum determination 
while assuring that there are no asymmetries or long tails, which further re- 
duces the danger of border effects and biases associated with a rapidly chang- 
ing spectral index. Again, it can be seen from Figure 5.b that the r opt method 
produce, at any given input energy very approximately normal distributions. 



4-2 Bias in the reconstructed energies 

Figure 6 shows the relative reconstruction error as a function of the injected 
energy for both reconstruction techniques. In the case of events without sat- 
urated detectors (Figure 6.b), there is no appreciable bias using r opt , while 
using ro there is an energy dependent bias which, at larger energies, can reach 
~ 10%. The difference is much more significant in the case of events with 
saturated detectors (Figure 6. a): the r opt approach produces almost negligible 
bias in the whole energy range while the reconstructed energy is largely un- 
derestimated using r . The later underlines the fact that r opt is very different 
for both populations of events. Analogous results are obtained for any array 
grid size. 

The energy reconstructions of events with and without saturated detectors 
are shown separately for both techniques in the scatter plots of Figure 7. As 
commented before, better reconstruction is achieved using r opt for both classes 
of events, while a more significant difference appears for events with saturated 
detectors. 
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Fig. 6. Bias in the reconstruction methods for events with (a) and without (b) 
saturated detectors. 
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Fig. 7. Reconstructed energy vs. real energy. Top: using S(ro). Down: using S(r op t). 
Left: events without saturation. Rigth: events with saturation. 

In fact, a main advantage of our proposed method is the treatment of events 
with saturated detectors, which is shown for r and r opt in Figures 7.b and 
7.d respectively. Using a fix characteristic value tq, events with saturation 
are poorly reconstructed, specially at lower energies. Essentially, the main 
problem is that these events have very few triggered stations and almost at 
the same distance from core. Consequently, their reconstruction accuracy is 
quite poor - this is particularly true for the fit to the LDF. In practice, using 
the r approach, these events probably would not pass the usual quality cuts 
and would be discarded, or would be reconstructed with a specific procedure. 
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Fig. 8. Fraction of events with saturated detectors as a function of energy for the 
different array spacings considered, (a) 9 = 30°. (b) 6 = 60°. 

Nevertheless, the r opt approach makes it possible to infer without almost any 
bias the energy of all the events, with an accuracy comparable to that attained 
for events without saturation. 

The advantage of a homogeneous treatment for all classes of events is further 
stressed by the fact that events with saturated detectors are in general dom- 
inant for most of the operational range of an experimental array, regardless 
of the detector separation (see Figure 8). For example, considering the 1 km 
separation array used for these studies, the number of triggered detectors in 
an event varies from 5 to 60 for showers from 10 17,5 to 10 19 ' 5 eV. Considering 
an incoming event of E ~ 10 19 ' 5 eV and a zenith angle of 9 ~ 30 degrees, any 
detector located at < 550 m from the shower axis would very likely be satu- 
rated. Thus, independently of the position of the core inside the elementary 
cell, almost 100% of the events have at least one saturated detector. At still 
higher energies, even 2 or 3 detectors could be saturated. Furthermore, for 
the same spacing, 50% of the events will have at least one saturated detector 
above E ~ 6 EeV (see Figure 8). 



4-3 Reconstruction of a rapidly changing spectrum 

In the previous section we demonstrated that the energy error distribution 
functions produced by the r method are wider, more skewed and have more 
extended tails than those produce by r opt . In principle, depending on the 
magnitude of these differences, they could affect the determination of spectral 
features, specially if the spectral index is varying rapidly over a narrow energy 
interval such as, for example, the ankle region and beyond. 

In order to assess the potential effects of using either technique for the re- 
construction of a structured spectrum with rapid changes as a function of 
energy, we use the following semi-analytical example. An idealized sectionally 
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continuous spectrum is assumed, that resembles the main spectral structures 
above 10 18 eV: the ankle, the GZK flux suppression [21,22] and a smooth low 
energy cut-off reflecting the discreteness of the surface array. The two latter 
suppressions in the spectrum are represented by functions of type tanh() of 
the input energy, while the remaining structures are represented by power laws 
separated by abrupt discontinuities in the first derivative. 

In order to reproduce analytically the energy error distribution functions given 
in Figure 5 as a continuous function of energy, we fit our simulation results 
with an Asymmetric Generalized Gaussian function (AGG): 



AGG 



(y) 



C7a 

k r(i/c) 
.r(i/c) 



ex p{-7/ c [-(z/ - /^)] c } if V< V 

exp{-7"[(y - ^)] c } if y>V 



where, 



<ti + a r \T(l/c) J a x \T(l/c) J a r \T(l/c) J 

and af and af are the variances of the left and right sides respectively of the 
probability density function and T(x) is the Gamma function. If af = af AGG 
is symmetric. Furthermore, if af = af and c = 2, AGG reduces to the regular 
Gaussian distribution function and, for c = 1, it represents the Laplacian 
distribution. 

The error functions determined previously in Section 4.1 have been fitted using 
the AGG function for the both techniques: tq = 600 m and the shower-specific 
r opt . In the latter case the fit reduces very nearly to a Gaussian distribution 
function while, for 7*0, the best simultaneous fit to the right and left a^s and 095 
C.L. shown in Figure 5, is obtained for c > 2.05. In this way we can reproduce 
the asymmetries present on the error distribution functions and analytically 
map real energies onto reconstructed ones over the whole energy range of the 
spectrum. 

The results are shown in Figure 9. It can be seen that, if both events with and 
without saturation are lumped together, the large wings associated with the 
error distribution functions of the r approach significantly distort the spectral 
features. In this particular example, the ankle, is widen and shifted, while the 
bump and GZK suppression are shifted upwards and much more pronounced. 
The r opt approach, on the other hand, fits very tightly the original spectrum 
with the exception of very low energies, near the full efficiency edge, due to 
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Fig. 9. Input model spectrum (black/thin line) and the reconstructed spectra using 
p(ro) (red) and p(r opt ) (blue/dashed line) as energy estimators. The right axis shows 
the fraction of isotropic events between 9 = 0° and 45° with saturation as a function 
of energy (thick dotted line) for the same array with 1000 m separation. 

border effects. The r$ approach can also give an equivalent fit, although noisier, 
if only events without saturation are used. However, the decrease in statistics 
by neglecting events with saturated detectors (cf., the fraction of events with 
saturation - Figure 9, right vertical axis) is so drastic that the reconstructed 
spectrum is only limited to a much shorter energy interval well before the 
GZK suppression. 

In order to understand the relative magnitude of these effects, one must note 
that at the AGASA experiment [18], for example, the systematic uncertainty 
in energy determination is around 18%. Three different sources of uncertainties 
combine to give this value. The first one is related with the detector, mainly 
its linearity (±7%) and response (±5%). Second, the uncertainties coming 
from the lack of knowledge of the LDF (±7%), the attenuation curve used 
(±5%), the shower front structure and delayed particles (which contribute 
±5% each). Finally the relation of p(ro) with energy (due to the hadronic 
interaction model supposed, simulation codes, chemical composition etc.), in- 
troduce an uncertainty of ±12%. In addition, they find an underestimation 
of 10% in reconstructed energies due to energy calibration with p(r ), which 
is compensated by the overestimation due to the shower front structure and 
delayed particles (5% each one). We proposed that the uncertainties related 
to the LDF and p(r ) determination could be reduced by using an r opt deter- 
mined on a shower to shower basis. However, while this may be a significant 
improvement, the other uncertainties would still dominate. 

The Auger Observatory, a hybrid detector, reports that the largest uncertain- 
ties [12] come from the fluorescence yield (±14%), the absolute calibration of 
FD (±10%), the FD reconstruction method (±10%). Systematic uncertainties 
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from atmospheric aerosols, the dependence of the fluorescence spectrum on 
temperature and on humidity are each at the 5% level. These uncertainties 
added in quadrature give a total uncertainty of 22% in fluorescence energy de- 
termination. Therefore, in addition to the fact that the method proposed here 
does not affect directly hybrid energy reconstruction because of the improved 
accuracy in the determination of the core position, the total uncertainty in 
the spectrum determination for hybrid experiments is widely dominated by 
FD uncertainties. 



5 Conclusions 

The primary CR energy is generally estimated in surface arrays by interpolat- 
ing the lateral distribution function of particles in the shower front at ground 
level at a fixed optimum distance tq from the shower core. This parameter is 
assumed to be predominantly dependent on the detector separation distance 
for a given layout geometry and, therefore, is considered as a constant for a 
given array. 

In this work we propose an algorithm to evaluate an equivalent, but shower- 
to-shower optimal distance, which we call r opt . We have performed a thorough 
analysis of the dependence of r opt on energy and zenith angle, and demon- 
strate that, contrary to reference [14], these are not negligible factors. In fact, 
not taking into account an event-specific r opt , produce wider error distribu- 
tion functions that can even affect the reconstruction of a highly structured, 
rapidly varying spectrum. The shower-to-shower r opt approach, on the other 
hand, is an unbiassed estimator of the CR primary energy, which produce 
also narrower, symmetric, almost Gaussian error distribution functions for en- 
ergy reconstruction. Those properties of r opt can additionally lead to much 
more reliable spectral reconstruction. The differences emerging from the two 
procedures, r vs. r opt , when applied to spectral reconstruction may have as- 
trophysical implications, specially in the coming era of improved precision. 

An important aspect of the r opt approach is that it has the additional ad- 
vantage of allowing the same unified treatment for events with and without 
saturated detectors; something that, in the ro approach is generally not possi- 
ble, requiring either the selection of events through quality cuts, or the sepa- 
rate reconstruction with different techniques of the two types of events. Since 
the fraction of events presenting saturation is a rapidly increasing function of 
energy, the later greatly reduces the effective energy range for spectral recon- 
struction in almost all practical situations. 

For practical application to real experiments, however, a proper calibration 
curve should be deduced specifically for r opt , which would further optimize it 
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as an energy estimator. 
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